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Abstract 

Reaction networks are systems in which the populations of a finite number of species evolve through predefined 
interactions. Such networks are found as modeling tools in many biological disciplines such as biochemistry, ecology, 
epidemiology, immunology, systems biology and synthetic biology. It is now well-established that, for small population 
sizes, stochastic models for biochemical reaction networks are necessary to capture randomness in the interactions. The 
tools for analyzing such models, however, still lag far behind their deterministic counterparts. In this paper, we bridge this 
gap by developing a constructive framework for examining the long-term behavior and stability properties of the reaction 
dynamics in a stochastic setting. In particular, we address the problems of determining ergodicity of the reaction dynamics, 
which is analogous to having a globally attracting fixed point for deterministic dynamics. We also examine when the 
statistical moments of the underlying process remain bounded with time and when they converge to their steady state 
values. The framework we develop relies on a blend of ideas from probability theory, linear algebra and optimization theory. 
We demonstrate that the stability properties of a wide class of biological networks can be assessed from our sufficient 
theoretical conditions that can be recast as efficient and scalable linear programs, well-known for their tractability. It is 
notably shown that the computational complexity is often linear in the number of species. We illustrate the validity, the 
efficiency and the wide applicability of our results on several reaction networks arising in biochemistry, systems biology, 
epidemiology and ecology. The biological implications of the results as well as an example of a non-ergodic biological 
network are also discussed. 
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Introduction 

Reaction networks are used as modeling tools in many areas of 
science. Examples include chemical reaction networks [1], cell 
signalling networks [2], gene expression networks [3], metabolic 
networks [4], pharmacological networks [5], epidemiological 
networks [6] and ecological networks [7]. Traditionally, reaction 
networks are mathematically analyzed by representing the 
dynamics as a set of ordinary differential equations. Such a 
deterministic model is reasonably accurate when the number of 
network participants is large. However, when this is not the case, 
the discreteness in the interactions becomes important and the 
dynamics inherently noisy. This random component of the dynamics 
cannot be ignored as it can strongly influence the system's 
behavior [8-10]. To understand the effects of this randomness, 
stochastic models are needed, and the most common approach is 
to model the reaction dynamics as a continuous-time Markov 
process, whose states denote the current population size. Many 
recent works have employed such stochastic models to study the 
impact of noise [11-14]. 



In stochastic models, the underlying Markov process (X(/)) (20 

is a pure-jump process whose state space S contains all the 
population size vectors that are reachable by the random dynamics. 
The probability distribution of (X(t)),> 0 evolves according to a 
system of linear ordinary differential equations (ODEs), known as 
the Chemical Master Equation (CME) or Forward Kolmogorov 
Equation [15]. The dimension of the system of ODEs is equal to 
the number of elements in the state space S, with each element 
representing a possible combination of reacting species abundanc- 
es. When S is finite and small in size, the CME can be solved 
analytically since it is simply a small and finite system of linear 
differential equations. However, for infinite state-spaces an exact 
solution to the CME is difficult to obtain except in some special 
cases [16,17]. Beyond these special cases, current methods often 
rely on truncating the infinite state-space to obtain finite 
approximations of the CME [18], and then resorting to efficient 
numerical methods for their solutions. Such methods include 
Expokit [19], which is based on Krylov Subspace Identification, or 
the backward Euler method proposed in [20], among others. Such 
an approach works well only for relatively small systems, as the 
curse-of-dimensionality renders the numerical solution of the 
truncated master equation of larger systems prohibitive. Never- 
theless, recent methods based on Tensor Train (TT) and 
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Author Summary 

In many biological disciplines, computational modeling of 
interaction networks is the key for understanding biolog- 
ical phenomena. Such networks are traditionally studied 
using deterministic models. However, it has been recently 
recognized that when the populations are small in size, the 
inherent random effects become significant and to 
incorporate them, a stochastic modeling paradigm is 
necessary. Hence, stochastic models of reaction networks 
have been broadly adopted and extensively used. Such 
models, for instance, form a cornerstone for studying 
heterogeneity in clonal cell populations. In biological 
applications, one is often interested in knowing the long- 
term behavior and stability properties of reaction networks 
even with incomplete knowledge of the model parame- 
ters. However for stochastic models, no analytical tools are 
known for this purpose, forcing many researchers to use a 
simulation-based approach, which is highly unsatisfactory. 
To address this issue, we develop a theoretical and 
computational framework for determining the long-term 
behavior and stability properties for stochastic reaction 
networks. Our approach is based on a mixture of ideas 
from probability theory, linear algebra and optimization 
theory. We illustrate the broad applicability of our results 
by considering examples from various biological areas. The 
biological implications of our results are discussed as well. 

Quantized Tensor Train (QTT) representations [21,22] show that 
for CME problems that admit bounded TT ranks, storage costs 
and computational complexity that grow linearly in the number of 
species may be achieved. These and other methods for the 
numerical solutions of the CME remain active topics of research. 

When S is infinite or very large in size, the most common 
approach for approximating the solutions of a CME is by 
simulating a large number of trajectories of the underlying Markov 
process (X(t)) t > 0 , and using the sample values of X(t) to estimate 
the distribution at time t. Such simulations are performed using 
Monte Carlo procedures such as Gillespie's stochastic simulation 
algorithm (SSA) or its variants [23-25]. Since the simulation time of 
SSA depends linearly on the number of reactions that occur 
during the simulation time period, these procedures can be 
cumbersome for large networks. It is well-known that the 
stochastic effects caused by the random timing of reactions 
become less important when the population size is large. The 
dynamical law of large numbers proved by Kurtz [26] shows that 
under an appropriate scaling relationship between the population 
size, reaction rates and the system size, the stochastic model of a 
reaction network converges to the deterministic model, as the 
system size goes to infinity. Under this scaling relationship, one can 
also approximate the stochastic dynamics with certain stochastic 
differential equations (SDEs) that are easier to simulate and 
analyze [27,28]. However, these SDE approximations can only 
work when the population sizes of all the species in the reaction 
network are large, which is often not the case. For a detailed 
survey on the topic of estimating the solution of a CME, we refer 
the readers to the paper [29] which contains an exhaustive list of 
methods for this purpose. 

In many biological applications, one in interested in analyzing 
the long-term behavior or stability properties of a reaction 
network. This is fairly straightforward for deterministic models 
because many tools from the theory of ordinary differential 
equations can be used for this analysis [30] . However, the stability 
properties of stochastic models for reaction networks are difficult 
to verify for the following reasons. Let us consider a stochastic 



reaction network whose dynamics is represented by the Markov 
process (X(?)),> 0 with state space S. The evolution of the 
distributions of this Markov process is given by (/?(0)(>o which is 
the solution of the CME corresponding to the reaction network. 
Heuristically, we regard the stochastic dynamics to be stable when 
the family of distributions (p(0)<>o i- s "well-behaved" with time. In 
this paper, we consider several notions of "well-behaved" 
dynamics. The strongest of these notions is the concept of ergodicity 
[3 1] which means that there exists a unique stationary distribution 
71 for the Markovian dynamics, such that p(t)—*n as /— >oo, 
irrespective of the initial distribution p(0). This is analogous to 
having a globally attracting fixed point in the deterministic setting. 
If S is finite, the process is ergodic if and only if it is irreducible, in 
the sense that all the states in S are reachable from each other. It is 
hence enough to check irreducibility of the process using, for 
example, matrix methods [32,33]. Contrary to this situation, our 
main interest in this paper is in analyzing the stability properties of 
stochastic reaction networks with an infinite state space S. Note 
that in such cases, irreducibility no longer implies ergodicity, since 
the trajectories of the Markov process may blow up with time (see 
the carcinogenesis example in the discussion section). In this 
regard, ergodicity cannot be considered as a generic property of 
reaction networks with infinite state-spaces since both ergodic and 
non-ergodic processes can be found in nature. Assuming ergodicity 
without verifying it beforehand seems to be therefore unreasonable 
from both theoretical and practical perspectives. The direct 
verification of stability properties like ergodicity is generally not 
possible as the CME cannot be explicitly solved, except in some 
restrictive cases [16,17]. The common approach of using Monte 
Carlo simulations for estimating the solutions of a CME is 
inadequate for assessing the long-term behavior and stability 
properties of a stochastic reaction network, because one can only 
simulate finitely many trajectories and those too for a finite 
amount of time. Some methods for analyzing stability properties 
without the need for simulations exist, but they either work for 
specific networks [17,34], very special classes of networks such as 
zero-deficiency networks [35], or assume system size approxima- 
tions where the stochastic dynamics is represented by an SDE 
[36,37]. Such system size approximations do not hold when some 
species are present in low copy numbers, and even if they hold, the 
approximation error generally blows up with time [28]. Hence the 
stochastic dynamics and the corresponding SDE may have 
completely different long-term behaviors. Our aim, in this paper, 
is to develop a theoretical and computational framework for 
analyzing the long-term behavior and stability properties of 
stochastic models for reaction networks that do not rely on 
computationally expensive Monte Carlo simulations or on system 
size approximations of the stochastic dynamics. A similar goal is 
also achieved in the works [38,39] where results on stability and 
moments bounds are also obtained. The approach proposed in 
[40] is built upon a Foster-Lyapunov criterion [31] and a 
quadratic Foster-Lyapunov function in order to estimate the 
location of the stationary distribution. In the same, yet different, 
spirit, the proposed approach also relies on a Foster-Lyapunov 
condition but using a linear Foster-Lyapunov function that allows 
us to establish ergodicity, moment bounds, moment convergence 
and the existence of attractive sets for moments. While the 
approach in [40] is fully computational, the one we propose is also 
theoretical and allows us to conclude on structural properties of 
classes of networks such as structural ergodicity, structural houndedness of 
moments and structural convergence of moments. Our approach relies on a 
mixture of simple ideas from stochastic analysis, linear algebra, 
polynomial analysis and optimization. Even though our conditions 
are only sufficient, we demonstrate their broad applicability by 
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successfully establishing stability properties of several reaction 
networks taken from the literature. 

We mentioned before that the stochastic and the deterministic 
models of a reaction network are connected through the 
dynamical law of large numbers [26]. It might be tempting to 
think that the stability properties of a stochastic model can be 
assessed by studying the stability properties of the corresponding 
deterministic model. However in general, the stochastic and 
deterministic models can have very different stability properties. 
This is because a deterministic model cannot capture noise 
induced effects which may have a significant impact on the long- 
term behavior of a system. For example, in the synthetic Toggle 
Switch by Gardener [41], the deterministic model exhibits 
bistability and hence starting from different initial values, the 
system can converge to two different steady states. On the other 
hand, the corresponding stochastic model is ergodic (see network 
(35)) and hence the solution of the CME converges to the same 
stationary distribution irrespective of the initial distribution. A 
similar phenomenon occurs with the repressilator (see [42] and 
network (36)), where the stochastic model is ergodic while the 
deterministic model exhibits oscillations. On the other hand, it is 
also possible to find networks for which the deterministic model 
has a locally asymptotically stable equilibrium point, implying 
that whenever the initial condition is contained within its region 
of attraction, the trajectories converge to it. If the initial condition 
lies outside this region of attraction, then the trajectories of such a 
network become unbounded with time. In the stochastic setting, 
the randomness causes each trajectory to leave the region of 
attraction in finite time, and then become unbounded suggesting 
that there is no stationary distribution for the dynamics (see 
network (22) and Figure 1). This lack of stationary distribution is 
because the stochastic dynamics can jump potential wells from 
one macroscopic fixed point which is stable to another fixed point 
which is unstable [43]. A more striking example of divergent 
deterministic and stochastic behaviors is given by network (26) 
(see also Figure 2). While the deterministic model has a unique 
globally stable fixed point, the stochastic model is non-ergodic 
and all the moments grow unboundedly with time. In this 
example it is impossible to predict the stochastic behavior from 




Time 



Figure 1. Trajectory of the state of the deterministic system 
(23) with initial condition Ko=0 (top); Sample path of the 
Markov process describing the network (22) with initial 
condition .\o = 0 (bottom). Whereas the trajectory of the state of 
the deterministic model converges to a stationary value, the trajectory 
of the state of the stochastic model goes unbounded. 
doi:10.1371/journal.pcbi.1003669.g001 
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Figure 2. Comparison of the trajectories of the deterministic 
and stochastic (first-order moments) models of the reaction 
network (22) with initial condition ;m(0) = 0, ;c 2 (0) = 2, JTi(0) = 0 
and \';{0)-2 for the deterministic (top) and stochastic dynam- 
ics (bottom), respectively. We can see that while the deterministic 
trajectories converge to their equilibrium point, the first-order moments 
grow without bound. 
doi:1 0.1 371/journal.pcbi.l 003669.g002 

the deterministic model. The above examples illustrate that the 
stability properties of the stochastic dynamics can, in general, not 
be assessed from the stability properties of the deterministic 
dynamics. 

Our results can help in understanding the stability properties of 
the moments of a Markov process (^(0)f>o representing a 
reaction network. In particular, we present a method to check if 
these moments remain bounded with time and if they converge to 
their steady state values as time goes to infinity. Such results can 
help in verifying the suitability of a model for a given system and 
in designing biological controllers that drive the moments to 
specific steady state values. We provide easily computable bounds 
for the moments that hold uniformly in time. We also determine 
bounds for the steady state moment values, which can help in 
understanding the properties of the steady state distribution, even 
if this distribution is not explicitly computable. In many biological 
applications, it is of great interest to explicitly compute the first 
few moments of the process (X(f)),>o without solving the 
corresponding CME. One can easily express the dynamics of 
these moments as a system of ordinary differential equations, but 
generally this system is not closed when the network has 
nonlinear interactions. Many moment closure methods that 
suggest schemes to close these equations to obtain approxima- 
tions for the moments have been proposed (see e.g. [44,45] and 
references therein). The results obtained in this paper can be used 
to ascertain the correctness of a given moment closure method for 
a specific network (see the example based on the network (29)). 
Furthermore, several moment closure methods are developed 
under an implicit assumption that the moment-generating function 
corresponding to the solution of the CME exists for all times. One 
of our results provides a way to easily check that this assumption 
is indeed valid. 

Reaction networks 

Let us now formally describe reaction networks. Motivated by 
the literature on chemical kinetics, we refer to the network 
participants as molecules which may belong to one of d species 
Si,...,Sd- There are K reactions in the network and for any 
k= 1, . . . ,K, the stoichiometric vector Ck = (Ck,l >■■ ■ Zk,d) denotes the 
change in the number of molecules in each of the species due to 
the fc-th reaction. 
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Deterministic models 

Consider the deterministic model for the reaction network 
described above. In this setting, the state of the system is described 
by a vector of concentrations of the d species which we denote by 
keR^q. The concentration of a species is simply its molecular count 

divided by the system volume. Let kk(K) be the flux associated with 
the k-th reaction (see [8]). To ensure positivity of the system, we 
require that 1a-(k) = 0 whenever K,- = 0 and Ck,i < 0. If the initial state 
is Ko, then the evolution of concentrations is given by (<l> Ko (t))t>o 
which satisfies the Reaction Rate Equations (RRE) of the form 

d(j> K (t) * . 

— fr- = E-W^WX* with ^ 0 (0) = Ko . (1) 

k=l 

We are interested in the long-term behavior and stability of our 
reaction dynamics. More precisely, we would like to check if the 
following conditions are satisfied. 

DC1 For any 7<o, there is a compact set AC(ko) such that 
^ KO WeiC( Ko ) for all f>0. 

DC2 There exists a compact set /Co such that for any Kg, we 
have (j> Kf) (t)elCo for large values off. 

DC3 There is a /c e9 such that for any Kq we have (j> Ko (t)— ►Keq as 

/->00. 

The first condition, DC1, says that for any Ko, the entire 
trajectory (<^ Ko (0)f>o stays within some compact set. We would 
expect this to be true for most realistic systems. Hence a violation of 
this property may suggest a flaw in the deterministic model. The 
second condition, DC2, says that there is an attractor set for the 
dynamics, where all the trajectories eventually lie, irrespective of 
their starting point. The last condition, DC3, says that there is a 
globally attracting fixed point for the deterministic model. Using 
techniques from the theory of dynamical systems [30,46] , one can 
verify these conditions, without the need of simulating the 
deterministic model. There is also a general theory to check 
condition DC3 for reaction networks satisfying mass-action 
kinetics (see [47-50]). Broadly speaking, these three conditions 
present different ways of saying that the reaction dynamics is 
"well-behaved". Our goal in this paper is to develop a theoretical 
and computational framework for verifying conditions similar to 
DCl, DC2 and DC3 for stochastic models of reaction networks. 

Stochastic models 

Consider the stochastic model corresponding to the reaction 
network described above. In this setting, the firing of reactions are 
discrete events and the state of the system refers to the vector of 
molecular counts of the d species. When the state is X, the A>th 
reaction fires after a random time which is exponentially 
distributed with rate Ak{x). The functions X\,...,)-k are known 
as the propensity functions in the literature. To ensure positivity of 
the system, we require that if x + Ck & > tnen ^k(x) = 0, where No 
is the set of non-negative integers. The dynamics can be 
represented by the Markov process (^Cv 0 (0)/>o where Xq is the 
initial state. Note that if X Xo (t) = (Xi(t), . . . ,X d (t)), then X^t) is 
the number of molecules of Si at time t. 

It is important to select a suitable state space 5 for the Markov 
process representing the reaction dynamics. We choose 5 to be a 
non-empty subset of satisfying the following properties: 

(A) If xeS and lk(x) > 0 for some k=l,...,K, then x + Cjt&S. 

(B) There is no proper subset Si a S satisfying part (A). 

Observe that part (A) ensures that if XqgS then X Xo ( t)eS for all 
/ > 0 and hence 5 can be taken to be the state space of all the 



Markov processes describing the stochastic reaction network with 
an initial state Xq in 5. Part (B) implies that the reaction dynamics 
cannot be contained in a proper subset of S. The role of this 
assumption will become clear in the next section, when we discuss 
the issue of state space irreducibility. Note that in certain cases, such 
as the pure-birth network 0^Sj, a suitable state space satisfying 
the above criteria cannot be found. There also exist cases where 
the above criteria restricts the choice of state space. For example, 
for the pure-death network Si — *-0, the only possible choice for state 
space is 5 = {0}. Finally we remark that if the reactions in a 
network satisfy a conservation relation then the state space must be 
chosen with an initial condition in mind. For example, for the 
network Si^±S2, the sum of molecular counts of Si and S2 is 
preserved by the reactions. Hence if we wish to study the stochastic 
dynamics with the initial sum as n, then the correct choice for state 
space is 5 = {(xi,X2)ef^l : X\ + x 2 = n}. 

Let V(S) denote the space of probability distributions over S, 
endowed with the weak topology which is metrized by the 
Prohorov metric (see [51]). For any x,yeS let p x (t,y) denote the 
following probability 

p x (t,y) = P(X x (t)=y). (2) 

Defining p x (t)(A) = v ^a Px(t,y), for any A a S, we can view p x (t) 
as an element in V(S). In iact,p x (t) is the distribution at time t of the 
Markov process (JG(0)«>o- The dynamics of p x (t) is given by the 
Chemical Master Equation (CME) which has the following form: 

d ^f 1 = J2(pAt,y-C l M k (y-Ck)-pAt,y)hLv)), (3) 
m k =i 

where p x (0,y)=l if x = y and p x (0,y) = 0 for all y^=x. Theoreti- 
cally, one can find p x (t,y) for any />0 and yeS, by solving this 
system. However this system consists of as many equations as the 
number of elements in S. Hence an explicit solution is only possible 
when S is finite, which only happens in very restrictive cases where 
all the reactions preserve some conservation relation. Typically, S is 
infinite and solving this system analytically or even numerically is 
nearly impossible, except in some restrictive cases (see [16]). From 
now on, we assume that 5 is infinite. 

The above discussion shows that at the level of distributions, we 
can view the stochastic dynamics (X X(l (t)) t > 0 as the deterministic 
dynamics (p Xo (t)),>o, which satisfies the CME. However, the major 
difficulty in analyzing this deterministic dynamics is that it occurs over 
an infinite dimensional space V{S). Nevertheless we can recast the 
conditions DCl, DC2 and DC3 in the stochastic setting as below. 

SCI For any Xq, there is a compact set K,(xq)<^V(S) such that 
p Xo (t)erC(x 0 ) for all />0. 

SC2 There exists a compact set K,q<^V{S) such that for any 
X(,eS we have p Xo (t)eK,o for large values of /. 

SC3 There is a neV(S) such that for any .To we have p Xo (t)—>n 
as t— >cc. 

Each of the above conditions give an important insight about 
the long-term behavior and stability of the stochastic dynamics. 
The first condition, SCI, says that for every ee(0,l) we can find a 
finite set A c <^S such that each p Xo (t) puts at least (1 —e) of its mass 
in A c . In other words, the probability that the state of the 
underlying Markov process at any time t is inside A c is greater 
than (1— e). We would expect this to be true for most realistic 
models. If condition SC2 holds then the evolution of distributions 
have a compact attractor set in V(S), where all the trajectories 
eventually lie irrespective of their starting point. This suggests that 
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in the long run, the family of processes {(X Xo (t)),> 0 : A "o£<S}, 
spend most of their time on the same set of states. The last 
condition SC3 says that the evolution of distributions have a 
globally attracting fixed point n. If this holds, then the Markov 
process representing the reaction dynamics is ergodic with n as the 
unique stationary distribution. For understanding the long-term 
behavior of a stochastic process, ergodicity is a desirable property 
to have. In the long-run, the proportion of time spent by any 
trajectory of an ergodic process, in any subset of the state space is 
equal to the stationary probability of that subset (see (12)). In other 
words, information about the stationary distribution can be 
obtained by observing just one trajectory for a sufficiently long 
time. Such a result can have important applications. For example, 
consider a culture with a large number of identical cells with each 
cell having the same reaction network. If we can show that this 
intracellular network is ergodic, then by observing the long-term 
reaction dynamics in a single cell, using for example, time-lapse 
microscopy, we can obtain statistical information about all the cells 
at stationarity. Conversely, ergodicity allows us to obtain the 
stationary distribution of a single-cell by observing the distribution 
over the population, using for example flow cytometry. 

In this paper we develop a general framework for checking 
conditions SCI, SC2 and SC3. However, the scope of our paper 
is broader than that. As mentioned in the introduction, we obtain 
easily computable bounds for the statistical moments of the 
underlying Markov process and investigate when these moments 
converge with time. We also present conditions for the distribution 
of the process to be light-tailed. 

Results 

Preliminaries 

In this section we discuss the main results of our paper. In 
particular, we explain how conditions SCI, SC2 and SC3 can be 
verified without having to simulate the trajectories of the Markov 
process representing the reaction dynamics. Intuitively, these 
conditions can only hold if the Markov process has a low 
probability of hitting states that have a very large size. In our case, 
the states are vectors in R and so we can measure their size by 
using any norm on R d . The central theme of this paper is to 
demonstrate that for many networks, long-term behavior can be 
easily analyzed by choosing the right norm for measuring the state 
sizes. This right norm has the form 

d 

IWI,= $>N, (4) 

1=1 

where v is a positive vector in U d satisfying the following condition. 

Condition 1 (Drift-Diffusivity Condition). For a positive 
vector veR"*, there exist positive constants C\,C2,c-$,C4 and a nonnegative 
constant C5 such that for all xeS 

K 

c\— C2<v,x> and (5a) 

k=\ 



K 

^2 4(-v)<v,U-> 2 <c 3 + f 4 <i',A-> + c 5 <v,x> 2 . (5b) 
k=l 

Here <y> denotes the standard inner product on U d . If we 
consider the process (|| JG 0 ML)r>o> then its dynamics can be seen 



to have two components drift and diffusion which have the form 

Sf=l 4W(i'.U> an d Y!k = \ 4W<v,4i> 2 respectively when 
X Xo (t) = x (see page 2 in the Supplementary Material SI). 
Condition 1 gives upper-bounds for the magnitude of these two 
components and hence we call it the drift-diffusivitj condition 
(abbreviated to Condition DD from now on; the abbreviations 
DDI and DD2 stand for the first and second inequality, 
respectively). Observe that when the process (||JGt 0 (0llv)r>0 goes 
above C\ jci then it experiences a negative drift, suggesting that it 
will move downwards. This fact will be crucial for our analysis. 

For now, we assume that a vector v satisfying Condition DD has 
been found. In later sections we demonstrate how v can be 
determined for a large class of networks by solving suitably 
constructed optimization problems. 

For any positive integer r, let >n r Xa (t) denote the r-th moment of 
\\X Xo {t)l defined by 

< o (0 = E(||Z, 0 (Oli;;) = £ \\y\\ r j> Xo (t,y). (6) 

yeS 

Similarly let *¥'(xo,t) denote the r-th moment of X XQ (t) at time t. 
Then *P'(.Yo,0 is a tensor of rank r whose entry at index 
. . . ,/, )e{l,2, . . . ,d}' is given by 

>l<;. Ax,,n Vr, ; ...r,,./ )>l: (/,r>. (7) 

where y = (yi, ■ ■ ■ ,yd) and p Xo (t) is the distribution of X Xfj (t). 
Suppose that for some positive constants r and C r (xo) we have 

sup m x (t)<C r (xo). (8) 

For any M>0, let Km be the compact (finite) set defined by 
Km = {xgS : ||x||,,<M} and let K° M denote its complement. 
Markov's inequality (see [52]) implies that for any e>0 we can 
choose M large enough to satisfy 

sup Pxo (t,K M )= sup p(h^ 0 (/)ii;;>m'-) < 
supe(k 0 (oii;)^< £ . 

Hence Prohorov's theorem (see Chapter 3 in [51]) ensures that 
condition SCI holds. Similarly we can prove that condition SC2 
will hold if for some r>0 there exists a constant C r such that 

lim sup m' x (t)<C r for all x 0 eS. (9) 

Relations (8) and (9) give uniform and asymptotic upper-bounds 
for m r Xo (t). Using these relations we can also obtain uniform and 

asymptotic upper-bounds for the entries of *P'(.Vo,?)- Such moment 
bound results have applications in queuing theory and control 
theory (see [53]). In Theorem 2 we show that under certain 
conditions, (8) and (9) hold and the upper-bounds can be easily 
computed. 

Instead of the r-th moment of the process (|| X Xo (t)\\ v ) t > 0 , one 
can ask if the exponential moment of this process is uniformly 
bounded from above. This will happen if for some y>0 we have 
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sup E(V I|Z; 

/>0 V 



> 0 Wllv 



sup 

'so ,, eS 



(10) 



If (10) holds, then the distribution p XQ (t) is light-tailed (a distribution 
is called light- tailed if its tails are majorized by an exponential 
decay) uniformly in /. This shows that all the cumulants of the 
distribution p Xfj {t) exist, which is an important result for the 
following reason. There is a considerable body of research 
dedicated to estimating the moments of the process (^x 0 ( z )) ) >o 
directly without computing the distribution functions p X(j (t). For 
any integer r > 0, one can easily write the differential equations for 
the dynamics of the first r moments. However when the reaction 
network has nonlinear interactions, this system of equations is not 
closed for any r. Various moment closure methods (see [54,55]) 
exist that specify ways to close these equations artificially and 
estimate the moments approximately. A popular moment closure 
method is the cumulant-neglect method which ignores the higher 
order cumulants of the distribution Px~(t) for all 7/>0. Of course 
this method is only valid when the higher order cumulants exist. 
This is guaranteed if (10) holds. In Theorem 3 we give conditions 
for verifying (10). 

We now come to the question of checking condition SC3 which 
says that the process (X Xo (t)) t > 0 is ergodic. This can only happen 
if the state space 5 is irreducible, which means that all the states are 
accessible from each other. Recall the definition otp x (t,y) from (2). 
Mathematically, we say that 5 is irreducible if for all x,yeS, we 
have p x (t\,y)>0 and p y (t2,x)>Q for some ?i,?2>0. In order to 
check the irreducibility of 5, one has to verify that there is no 
proper subset 5i <=5, such that once the process reaches a state in 
S\ , it stays in 5i forever. For reaction networks with mass-action 
kinetics, methods for checking irreducibility have recendy been 
reported in [56] and [57]. These methods can be extended to 
situations where the propensity functions are positive in the 
positive orthant. When the propensity functions vanish inside the 
positive orthant, the problem of checking irreducibility can 
become much more complicated, and to the best of our knowledge 
no methods exist in the literature for this purpose. 

We mentioned before that the vector v is chosen so that the 
process (|| JG 0 (OL)r>o has a negative drift at large values. 
Assuming irreducibility, this is sufficient to verify ergodicity of 
(X Xo (t)),>o (see Proposition 4). 

Suppose that condition SC3 is satisfied and the process 
(JG 0 (0)f>o is ergodic with stationary distribution n. For any 
positive integer r, let IT denote the r-th moment of the stationary 
distribution n. Then TV is a tensor of rank r defined in the same 
way as *¥ r (x<j,t) (see (7)), with p x „(t,y) replaced by n(y). Using 
Theorem 2 we can determine the values of r for which II' is finite 
(componentwise) and W ' (xo,t)—*Tl r as t— >oo (see Theorem 5). We 
can also identify functions / : S— >R for which 

lim Ef/(X, 0 (f)))=y70)7tOO<oo (11) 

yeS 

holds for any x$eS. If / is such a function, then the ergodic 
theorem for Markov processes (see [58]) says that 



lim - 

t^OJ t 



f(X XQ (s))ds = $3/60*00 almost surely, (12) 



for any x^eS. Lastly, we also obtain conditions to check if the 
stationary distribution n is light-tailed (see Theorem 6). 



General results 

In this section, we formally present the main results of our 
paper. Their proofs are given in the Supplementary Material SI. 

Moment bounds. Our first result establishes that for certain 
values of r, we can obtain uniform and asymptotic moment 
bounds for the r-th moment of the process (||ix o (0L)/>o- 

Theorem 2 Assume that Condition DD holds. Let r max he given by 



1+^ ifc 5 >0 
oo if cs= 0. 



(13) 



For any positive integer r, if r < r max then there exist positive constants 
C r (xo) and C, such that (8) and (9) hold. 

The values of the constants C,(xq) and C r can be explicidy 
computed using a recursive relationship (see the Supplementary 
Material SI). Note that if V = (vi, . . . ,vf), then for any 
y = (}>l, ■ ■ ■ ,yd)sS we have >',< ||v||,,/v, for any i. Hence for any 
i u ...,i r e{\,2,...,d} we have" % ir (x 0 ,t)<m' Xo (t)/ W j=l Vij 
Therefore using Theorem 2, we can obtain uniform and 
asymptotic moment bounds for the reaction dynamics 
(X Y(1 (/)), 20 (see the Supplementary Material SI). 

Observe that if C5 = 0 then r max = 00 . In this case, Theorem 2 
says that for each positive integer r and x$eS there exists a 
constant C r (xo) such that (8) holds. By showing that we have a 
C>0 such that C,(.Vo)<r!C for all positive integers r, we obtain 
our next result, which gives sufficient conditions to check (10). 

Theorem 3 (Uniform Light-Tailedness) Suppose that 
Condition DD holds with C5 = 0. Given an initial state XqeS there exists 
a V>0 such that 



sup E 

f>0 



sup\_ 



v Px 0 (t,y)<<x>. 



yeS 



Ergodicity and Moment Convergence. The next result 
verifies the ergodicity of a reaction network satisfying Condition 
DD. It follows from Theorem 7.1 in Meyn and Tweedie [59]. 

Proposition 4 (Ergodicity) Assume that the state space S of the 
Markov process (X Xo (t)) t> Q is irreducible and Condition DDI holds. Then 
this process is exponentially ergodic in the sense that there exists a unique 
distribution neV(S) along with constants B,c>0 such that for any x$eS 



sup Jt,A)- 



■n{A)\<Be- c 'for all t>0. 



This result says that as t— >oo, the distribution p Xlj (t) converges 
to 71 exponentially fast. Henceforth we assume that the process 
(X vo (?)) ( >o is ergodic with stationary distribution n. 

Let / : iS— >IR be a function such that for some positive integer 
r<(r m sx— 1), there exists a C>0 satisfying \f(x)\ < C(l + [|jc|Q 
for all xeS. Using Theorem 2 we can prove that for such a /, the 
relations (11) and (12) hold. As a consequence we obtain the 
following result about the convergence of moments with time. 

Theorem 5 (Moment Convergence) Assume that Condition 
DD holds. Let r be any positive integer satisfying r < (r mllx — 1 ). Then TV is 
finite (componentwise) and (xo,t)^>Tl' as t—>co. 

If f(x)= \\x\\ r r then Theorem 2 and (11) imply that for any 
positive integer r<(r max — 1) there exists a positive constant C r 
such that 
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^lbll>0)<c r . (14) 

yeS 

In particular, if C5 = 0 then r max = go and ( 1 4) holds for each r. By 
proving the existence of a constant C>0 such that C r <r\C r for 
all positive integers r we get our last result which shows that the 
stationary distribution is light-tailed. 

Theorem 6 (Light-Tailedness at stationarity) Suppose that 
Condition DD holds with C5 =0. Then there exists a y>0 such that 

e rMv n(y)<co. 

yeS 

The framework described above is very general and can be 
applied to any network that satisfies Condition DD. In what 
follows, we specialize the results for two wide classes of networks 
with mass-action kinetics, namely reaction networks with mono- 
molecular and bimolecular reactions. It will be, however, pointed 
out in the examples that the scope of our approach is much 
broader since more general propensities, such as those involving 
Hill functions or more general mass-action kinetics, can be 
considered. 

Methods 

Using the analytical tools developed in the previous sections, 
several general results can be stated for the class of unimolecular 
reaction networks and bimolecular reaction networks. In what 
follows, when we say that a moment is bounded, we mean that it is 
bounded uniformly in time (as in (8)). This can be established using 
Theorem 2 once Condition DD is verified. Furthermore, when we 
say that a moment is globally converging, we mean that it 
converges to its equilibrium value as time tends to infinity, 
irrespective of the initial state Xq. Once, Condition DD is verified, 
this can established using Theorem 5. 

The main aim of the section is to develop a theoretical and 
computational framework for checking Condition DD. 

Results for stochastic unimolecular reaction networks 

Let us then consider a unimolecular reaction network which 
involves d species that interact through K reaction channels of the 
form: 

0 - 1 Si, Si - 0, Si -i- £vfSj (15) 

./=! 

where i= I, . . . ,d, £e{\, . . . ,N,}, N,>0 and vfeN 0 - The reaction 
rates k' 0 , k® and kj are positive real numbers. In accordance with 
(3), the reactions are indexed from n = 1 to K, and corresponding 
propensities and stoichiometries are denoted by ). n {x) and f„, 
respectively. 

Motivations. The unimolecular case may seem quite restric- 
tive at first sight and not of particular practical interest. We 
demonstrate below that, on the contrary, the proposed results on 
unimolecular reaction networks complete existing ones and are, 
therefore, of practical and theoretical interests. Although some 
explicit solutions for the CME are indeed known for some particular 
unimolecular reactions [16], it is still unknown whether the CME 
admits an closed-form solution for all possible type of unimolecular 
reactions. Note that no simplification nor assumption is made on the 



problem in our work. Hence, we are dealing with the very general 
unimolecular case. 

The results developed of this section are useful in several ways. 
First of all, all types of unimolecular reactions can be handled with 
the proposed approach, making it more general than existing ones 
in this regard. Moreover, given a specific reaction network, the 
method allows one to establish whether a unique stationary 
distribution exists without solving the CME. This is particularly 
important since unimolecular networks may not be ergodic. In this 
case, the network can exhibit unstable behaviour which may 
suggest a flaw in the model if the considered real-world system 
exhibits stable trajectories. Moreover, in certain design applica- 
tions such as those in synthetic biology, it seems natural to design 
networks that have well-behaved dynamics. Checking ergodicity 
provides a convenient way to determine if the network dynamics is 
well-behaved. Note, furthermore, that it is, in general, difficult to 
infer ergodicity directly from the solution of the CME (when it is 
known) since proving the existence of a unique globally attractive 
stationary distribution amounts to checking convergence of the 
solution to the CME to the same distribution for all possible initial 
distributions, which are in infinite number in our setup. This fact is 
even more true when large networks are considered since the 
explicit form of the solution to the CME is, in this case, very 
intricate [16]. The proposed results allow one to circumvent this 
difficulty and demonstrate that ergodicity can be assessed by very 
simple means, i.e. using basic notions of linear algebra. The results 
can be furthermore used to assess the structural ergodicity of a 
reaction network, that is, the ergodicity of a network for any 
combination of the rate parameters, by very simple means. This 
very strong and practically relevant notion is extremely difficult, 
again, to check from the solution of the CME since it would 
require to check the convergence of the solution of the CME to the 
same stationary distribution for all initial conditions and all 
positive values of the rate parameters, a very cumbersome task, 
even for small networks. Finally, the results pertaining to 
unimolecular networks will also turn out to play an important 
role in the ergodicity analysis of bimolecular reaction networks. 

Theoretical results. Let us start with several theoretical 
results that characterize the long-term behavior of unimolecular 
networks of the form (15). 

Proposition 7 (Ergodicity of unimolecular networks) 
Let us consider the general unimolecular reaction network (15) and assume that 
the state-space of the underlying Markov process is irreducible. Let the matrices 
AeU dxd andbeU<l 0 , \\b\\^Q, be further defined as 

K 

A„(x)<v,L> = x T Av + b T v. (16) 

«=1 

Then, the following statements are equivalent: 

1 . The matrix A is Hurwitz-stable, i.e. all its eigenvalues lie in the open left 
half-plane. 

2. There exists a vector veR'^Q such that Av<0. 

Moreover, when one of the above statements holds, the Markov process 
describing the reaction network is exponentially ergodic and all the moments are 
bounded and globally converging. 

The above result shows that, for unimolecular networks, 
ergodicity and the existence of moment bounds can be directly 
inferred from the properties of the matrix A defined in (16). The 
second statement, which characterizes the Hurwitz-stability of A 
in an implicit way, will turn out to play a key role in the analysis of 
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unimolecular and bimolecular reaction networks since checking 
whether Av<0 for some v>0 is a linear programming problem. 

It is important to stress that, in the result above, if we simply 
demand that the moments be bounded and converging, then A 
may be allowed to have zero eigenvalues in certain cases. Note, 
however, that the moments will converge to values that may 
depend on the initial conditions. 

In the case that the structure of the network (the reactions and 
the stoichiometrics) is exactly known, but that the reaction rates 
are subject to uncertainties, the above theorem can be robustified to 
account for these uncertainties. To this aim, suppose that the 
matrix A depends on a vector <5e[ — 1,1]'' where f/eN is the 
number of distinct uncertain parameters. We write this matrix as 
A(S) and assume that there exists a matrix A + eM. lix ' t satisfying 
the following properties: 

1. A(S)<A + (in the componentwise sense) for all Se[— 1,1]'' 

2. There exists a S*e[— 1,1]'' such that A + =A(S*). 

Note that such a matrix A + may not exist, especially when 
some entries are not independent. However, when A + exists we 
have the following result. 

Proposition 8 (Robust ergodicity) Let us consider the general 
unimolecular reaction network (15) described by some uncertain matrices A(3) 
and b(S), \\b(5)\\^0. Assume further the matrix A(S) admits the upper- 
bound A + defined above and that the state-space of the underlying Markov 
process is irreducible for all uncertain parameter values (5e[— 1,1] . Tlien, the 
following statements are equivalent: 

1. The matrix A{5) is Hurwitz-stable for all 8e[— 1,1]''. 

2. The matrix A + is Hurwitz-stable. 

3. There exists a positive vector veR d such that A + v<0. 

Moreover, when one of the above statements holds, the Markov process 
describing the reaction network is robustly exponentially ergodic and all the 
moments are bounded and globally converging. 

Observe that checking the Hurwitz-stability property of each 
A(S) is equivalent to checking it for only A + . Hence we can 
conclude that, in this case, checking ergodicity of a family of networks is 
not more complicated than checking ergodicity of a single network. The case 
when the matrix A + is not defined is more complicated and is 
discussed in the supplementary material SI. 

Computational results. We now present several computa- 
tional results that accompany the theoretical results of the previous 
section. It is possible to extract many computational results from 
our general framework, but for simplicity we only address the 
problems of checking ergodicity and computing the first-order 
moment bounds. The asymptotic first-order moment bound, 
defined in Theorem 2, is given by C\=C\/c2. So the question 
arises: what is the smallest value for such a ratio? Or, in other 
words, what is the smallest attractive compact set for the first-order 
moment of <v,X(f))? Several numerical methods, solving exactly 
or approximately this problem, are discussed in the supplementary 
material SI. One of them is the following optimization problem 
which is fully equivalent to Proposition 7: 

Optimization problem 9 Let us consider the general unimolecular 
reaction network (15) and assume that the state-space of the underlying 
Markov process is irreducible. Assume further that the optimization problem 

max z ,,z s.t. z>0,v>£ 

17 

(zI + A)v<0 K ' 

is feasible with (z*,v*) as minimizer. Then, we have < b T v* / z* and 
Proposition 7 holds. 



A striking feature about the above optimization program is that 
the numbers of variables and constraints are given by d + 1 and 
2d+l, respectively. This means that the optimization problem 
scales linearly with respect to the number of species (d) in the 
network, and is independent of the number of reactions K. 
Therefore, from the point of view of this optimization problem, the 
size of a unimolecular network can be identified with the number 
of species, and not the number of reactions. The above 
optimization problem can be efficiendy solved using a bisection 
algorithm over z that is globally and geometrically converging to 
z* . Each iteration consists of solving a linear program, a class of 
optimization problems known to be very tractable, and for which 
numerous advanced solvers exist [60] . These properties, altogeth- 
er, make the overall approach highly scalable, which is necessary 
for dealing with very large networks. 

Results for stochastic bimolecular reaction networks 

Similar results are now presented for stochastic bimolecular reaction 
networks which, in addition to the unimolecular reactions (15), also 
involve bimolecular reactions of the form: 

Si + Sj ELivfS m , Si + Sj 1 0 ( 18 ) 

where ij = l,...,d, £e{\, . . . ,Ny}, Ng>0, and vfeN 0 . The 
reaction rates k\- and are positive real numbers. 

Theoretical results for bimolecular networks. When 

bimolecular reaction networks of the form (15)— (18) are considered, 
the left-hand side of condition (5a) can be expressed as 

K 

h(xKvX k ) = x r M(v)x + x J Av + b J v (19) 

;=1 

where M( v)eU d x d is symmetric, AeU dxd and beUt, 0 . Let 
S : = [ £i ... Cjj ] be the stoichiometry matrix of the bimolecular 
reaction network (1 5) — (1 8), and let S q be the restriction of S to 
bimolecular reactions, only. Further define a set 

J\fg : ={veR d : v>0, i^S^O}. 

When veNq, the quadratic term x T M(v)x in (19) vanishes, and 
equation (19) reduces to 

K 

Y,h(xKv,C k } = x T Av + b T v 
1=1 

which is exactly the same expression as in the case of unimolecular 
networks. This means that, with the additional constraint that 
veM q , all the results derived for unimolecular networks direcdy 
apply to bimolecular networks as well. This allows us to obtain the 
following result. 

Proposition 10 (Ergodicity of networks) Let us consider the 
bimolecular reaction network of the form (15)— (18) such that \\b\\=£Q in (19) 
and assume that the state-space of the underlying Markov process is irreducible. 
Assume further that the network admits a non-empty Afq. 

If there exists a vector vej\f q such that the inequality A v < 0 holds, then the 
stochastic bimolecular reaction network (15)— (18) is ergodic and all the 
moments are bounded and globally converging. 

It is important to mention that the existence of a non-empty 
set Nq is a prerequisite for utilizing the above result. Non- 
emptiness of J\fq is equivalent to the existence of a conservation 
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relation for all the bimolecular reactions, i.e. the value of (at 
least) a positive linear combination of the species populations 
remains unchanged when any of the bimolecular reactions fires. 
Note that this definition extends to more general mass-action 
kinetics as well. A necessary condition for the non-emptiness of M q 
is that S q is not full-row rank. This non-emptiness condition may 
seem restrictive at first sight, but it will be shown that several 
important reaction networks from the literature satisfy this 
condition. 

Whenever Af q is empty or there is no vej\f g such that Av<0 
holds, the next result can be used. 

Proposition 1 1 (Ergodicity of bimolecular networks) 

Let us consider the bimolecular reaction network of the form (15)— (18) such 
that ||/)||t^0 in (19) and assume that the state-space of the underlying 
Markov process is irreducible. Assume further that one of the following 
statements holds: 

1. There exists veU. d >0 such that Av<0 and M(v)<0 hold. 

2. There exists velR^g such that M(v) is negative definite. 

Then, the stochastic bimolecular reaction network (15)— (18) is ergodic and 
all the moments up to order (|_1 +2C2/C5J — 2) are bounded and globally 
converging. 

In the above result, the first statement can be checked using a 
linear program since the inequalities are componentwise. Check- 
ing the second statement, however, requires a semidefinite 
program, which is a more general convex program, that can be 
solved using solvers such as SeDuMi [61] and SDPT3 [62]. More 
details on the above result can be found in the supplementary 
material SI. 

Computational results for bimolecular networks. It is 

shown here that, once again, the theoretical results can be easily 
turned into linear programs that can be checked in a very efficient 
way. The following result is the numerical translation of 
Proposition 10. 

Optimization problem 12 Let us consider a bimolecular reaction 
network (15)— (18) and assume that the state-space of the underlying Markov 
process is irreducible. Assume further that Afq ^ 0 and that the optimization 
problem 

max Zj ,,z s.t. z>0,v>e 

(zI + A)v<0 (20) 

V T Sq=0. 

is feasible with (z*,v*) as minimizer. Then, we have < b T v* / z* and 
Proposition 10 holds. 

The computational complexity of this optimization problem 
scales linearly with the number of species and can therefore be 
solved for large networks. The following optimization problem is 
the computational counterpart of the first statement of Proposition 
11. 

Optimization problem 13 Let us consider a bimolecular reaction 
network of the form (15)— (18) and assume that the state-space of the 
underlying Markov process is irreducible. Assume further that the nonlinear 
optimization problem 

max zv z s.t. z>0,v>£ 

(zI + A)v<0 (21) 
M(v)<0. 



Long-Term Behavior of Stochastic Reaction Networks 

is feasible with (z*,v*) as minimizer. Then, we haveC\<b v*/z* and 
Proposition 11 holds. 

The above optimization problem does not scale as nicely as (20) 
since, in the worst case, the number of constraints related to M( v) 
is quadratic in the number of species. The problem, however, 
remains tractable due to the linear programming structure. 

Qualitative differences between deterministic and 
stochastic dynamics 

In this section we illustrate that stochastic and deterministic 
models of the same reaction network may exhibit very different 
qualitative behaviors. Therefore assessing ergodicity or the 
convergence of moments of a stochastic model from the stability 
properties of the corresponding deterministic model is, in general, 
incorrect. To support this claim, we consider two reaction 
networks. 

Jumping potential wells. Our first example shows that 
stochastic dynamics can jump potential wells and leave the stability 
regions of the deterministic dynamics, resulting in an unstable 
behavior. Consider the following reaction network: 




S 



S "i" 0 (22) 
S + S -i- 3S 

where 0 < a < /J. The deterministic dynamics for this network is 
given by 

k = f(K) : =k 2 -{i + P)k + aft (23) 

where keIR>o denotes the concentration of S. The fixed points for 
the dynamics are K_ = a and K + = /?, respectively. From the graph 
{(k/(/c))sIR>o x R : ;ce[R>o}, it is immediate that the fixed point 
K_ = a is locally asymptotically stable with the region of attraction 
as [0,/?) while the other fixed point K + =[i is unstable. 

We now consider the stochastic version of this network and let A 
be the generator of the corresponding Markov process. For the 
identity function f(x) = x we have 

Af(x) ^i^-L+p+^j x + *p. (24) 

The polynomial on the right-hand side has two positive roots that 
are 



x±=a+j8+i±y(a+0+^ -2a0. (25) 

This means that for all xeMrj satisfying ,v>l+.v + , we have 
f\f(x)>£, for some £>0, implying that the drift is positive. So if 
the state of the state of the network reaches a value that is greater 
than l+x + , then there is a possibility that the trajectories become 
unbounded with time. 

To demonstrate this, we pick a = 7/2 and /? = 21/2. In such a 
case, the largest root of the polynomial on the right-hand side of 
(24) is x+ = (29 + x/547)/2~26.194>j3. We can see that the 
region where the drift A/(x) is negative is actually larger than the 
region of attraction of the locally asymptotically stable fixed point 
for the deterministic dynamics. This is due to the fact that the 
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propensity function of the bimolecular reaction differs from 
whether we are in the deterministic or in the stochastic setting. 

Let us now set the initial condition Kq = 0 for the deterministic 
model and Xq = 0 for the stochastic one. Note that they both lie 
within the region of attraction of the fixed point of the 
deterministic dynamics and in the region of negative drift for the 
stochastic dynamics. We then perform 1000 SSA runs over 
100 seconds and stop the simulation when the propensity function 
x(x— l)/2 of the bimolecular reaction exceeds the value 
corresponding to 15000 molecules (approx. 1.12 xlO 8 ). At this 
rate value, the bimolecular reaction fires, on average, every 10~ 8 
seconds, leading to an explosion of the state of the system and to 
unbounded trajectories. Out of 1000 SSA runs, all were stopped 
before the end of the simulation time-period (100 seconds). This 
behavior strongly indicates that the system is not ergodic despite the 
fact that the deterministic model has a locally asymptotically stable 
fixed point. Figure 1 illustrates the above discussion. 

Globally stable deterministic dynamics does not imply 
moments stability. In the previous example, the stochastic and 
deterministic behaviors were different, but one can still understand 
stochastic instability through the deterministic model. The 
deterministic dynamics possesses a region in which the solutions 
explode and the randomness in the stochastic dynamics allows it to 
enter this region in finite time and grow unbounded thereafter. We 
now present an example which is more striking in the sense that 
the deterministic model cannot be used in any way to infer the 
instability of the stochastic model. In this example, the determin- 
istic dynamics has a unique fixed point which is exponentially 
stable, while the stochastic dynamics is not ergodic with all its 
moments growing unboundedly with time. 

Consider the reaction network given by 

0 1 Si 

0 -i S 2 (26) 
Si+S 2 -i 0. 

Let fce[R> 0 be the vector of concentrations. The state of the 
deterministic model evolves according to 



Finding an attractive compact set for the first-order 
moments 

The goal of this section is to compute a compact set that is 
attractive for the first-order moment of <v,X(/)> using the 
optimization problems (17) or (20). Due to the moment closure 
problem [54], analytical expressions for the steady-state values of 
the moments of bimolecular reaction networks are not available, 
and hence this is an important class of networks to analyze. 
Consider the following bimolecular reaction network 

0 - Si, Si - 0 
Si+Si *2 S 2 , S 2 % S, + S, (29) 
S 2 - 0. 

representing a dimerization process, i.e. Si dimerizes to S 2 . It is 
easily seen that this network is irreducible since any point in the 
state-space can be reached from any other point in a finite number 
of reactions having nonzero propensities. Choosing v in A/" ? , e.g. 
v T = [l 2], yields that c\=k and c\ = min{y 1 ,y 2 }, hence the 
network is exponentially ergodic, and all the moments are 
bounded and converging. On solving the optimization problem 
(20) with numerical values k=\, yi=y 2 = 0.2, k\2 = \ and 
fc21=0.1, we get that C\=c\/ct = S which coincides with the 
theoretical value A;/min{}'i,y 2 }. One can regard 
{(x\,X2)eK. 2 >0 : v t .y< Ci} to be an attractive compact set in 
which the first-order moments of <(v,X(/)> eventually lie. To 
validate this calculation, Monte-Carlo simulations were performed 
which yield 

lim E[<v,X(?)>] = 5.024 + 0.05, (30) 

t— *00 

showing the correctness of the attractive compact set. To further 
illustrate this result, several trajectories of E[Xi(?)] and E[X 2 (?)] for 
different initial conditions are plotted in Figure 3. 

We now discuss how the computation of an attractive compact 
set for the first-order moments can be used to assess whether a 
closure method leads to a result that is consistent with the 



Ki(0 = \-K\(t)K 2 (t) 
k 2 (t) = \-K X (t)K 2 (t). 

Assume that the initial conditions satisfy K 2 (0) — K\(Q) = a, for 
some aeR. Then we have the following result. 

Theorem 14 The unique equilibrium point of the dynamics (27) given 

by 

k\ = ^ \/a 2 + 4) and = ^ (W\/a 2 +4) . (28) 

is globally exponentially stable. 

In the stochastic setting, the picture is completely different as the 
next result indicates. 

Theorem 15 The Markov process corresponding to the stochastic model 
of network (26) is not ergodic and all its moments grow unboundedly with time. 
Moreover, if X\ (0) — X 2 (0) = a for some tx>0, we have that 
E{X\(t)~X 2 (t)]=(X for all t>0. 

To illustrate this result, we simulate the deterministic and the 
stochastic process (10000 SSA runs) for Ki(0) = 0, Ki(0) = a, 
X\(0) = 0, X 2 (0) = a and a = 2. The results are shown in Figure 2. 




Mi(t) 



Figure 3. Trajectories of the first order moments n 1 (i) = E[Xi(i)] 
and //,(/) -F,Vi(n of network (29) for different initial condi- 
tions (averaging is performed over 5000 cells). The trajectories 
converge to the unique steady-state value located inside the compact 
set (the surface below the dashed line), very close to the boundary. 
doi:1 0.1 371/journal.pcbi.l 003669.g003 
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stochastic dynamics. The idea is to check whether the closed 
system converges towards a value which lies within the compact 
set. Let us consider the reaction network (29) and close the first- 
order moments equations by neglecting the second order 
cumulant, i.e. neglecting the variance. By doing so, we get the 
model 



«i(0 = fc-yiM0-*i2/«i(0G'i(0-i)+2fc2iM2(0 

fa(t) = knM^i (0-1) -72^(0 



(31) 



where Jl\ and Jl 2 are the approximate first-order moments of the 
system. The unique positive equilibrium point for this model is 
given by 



Mi 

Pi 



i 

2k 



-Yv 



k\ 2 y 2 



2(y 2 +k 2 i) 



Yz+ kn 
«T(A«i-l) 



(32) 



where A = 



-Yl + 



ku7 2 Y,4kk u y 2 



y 2 + k 2 ij y 2 +k 2 i 
With the same parameter values as before, we find that 
^1 = 1 .6238 and /i^ = 1 .688 1 and therefore v T u* = 5 for 
v T = [ 1 2], showing that the state of the closed system converges 
to the boundary of the compact set. Note that SSA also predicts 
that the trajectories of the first-order moments converge to the 
boundary of this set. However the actual equilibrium values for the 
first-order moments of the stochastic dynamics are /4 — 1.1450 
and /( 2 — 1 .9350, which differ from the ones obtained with the 
closure method. This discrepancy is expected since the variance 
has been neglected. 

This example shows how attractive compact sets for the 
moments can be used as a test for the moment-closure methods 
by checking whether the closed system predicts trajectories that 
converge inside those compact sets. However, note that in the 
current state, these compact sets can only be used to obtain a lower 
bound on the closure-error whenever the trajectories of the closed 
dynamics converge to a point outside the compact set. In such a 
case, the lower bound on the closure-error e is simply given by the 
distance between the equilibrium point of the closed-system 



£> inf 112*- 



(33) 



where C is the attractive (convex) compact set and 2* is the 
equilibrium point of the closed dynamics. 

Feedback loop 

Let us consider the feedback loop network of Figure 4 
represented by the reaction network 



S. 

A '32 
03 — *■ 

Si ^ 



Si +s 2 , 
s 2 +s 2 , 



rx /(S3) « 

0 ^ Si 



(34) 



where Si is mRNA and S 2 is the corresponding protein. The 
dimer S3 acts back on the gene expression through an arbitrary 
bounded nonnegative function /(•). 
We have the following result: 



U 

AAAAA 



y 3 



-> 0 



Y 2 



"> 0 



-> 0 



|_> I J ^3> 



Figure 4. Feedback loop with arbitrary feedback rule. 

doi:10.1371/journal.pcbi.1003669.g004 



Result 16 For any positive values of the rate parameters and any bounded 
nonnegative function /(■), the feedback loop with dimerization (34) if ergodic 
and all the moments are bounded and globally converging. 

Stochastic switch 

Let us consider the stochastic switch of [63] described by the 
unimolecular stochastic reaction network 



(35) 



Above Sj and Sj represent mRNAs and proteins of gene i, 
respectively. The functions /i(-) and/iO) are arbitrary bounded 
nonnegative functions. We have the following result: 

Result 17 For any positive values of the rate parameters and any bounded 
nonnegative functions /i(') and f 2 (), the stochastic switch (35) is ergodic and 
all the moments are bounded and globally converging. 

Repressilator 

We consider here the stochastic repressilator of Figure 5 (see 
also [42]) involving N genes. The reaction network corresponding 
to this A^-gene repressilator is given by 



0 




s?, 


S? 


^ s? + s| 


0 


/ 2 <s|) 


si 


s» 


^ s?+s 2 


sj 




0- 
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0 


/i(^) 


sj 


0 




si 


0 


/ 3 <si) 


si 


0 




1 ! 


si 


k A 


s|+s? 




k A 


s 1 +s 2 

2 2 


»3 


h 


Og T" 03 


S 1 
a N 




Sjj + S^ 


sj 




0, f=l, 


sr 




0, i=l, 



(36) 



where f i (x) = a i + P t /(l+X n ), ot / ,j8 i ,«>0. Above, S? and S? 
are the mRNA and protein corresponding to gene i. We have 
the following result: 

Result 18 For any positive values of the rate parameters /r/ ,y, , i5, , tx, , yCf, 
and n, the stochastic N -gene repressilator (36) is ergodic and all the moments 
are bounded and globally converging. 



proteins A, R and C are depicted in Figure 7 where we can 
observe the expected oscillatory behavior. When averaging the 
populations of the proteins A, R and C over a population of 2000 
cells, we obtain the sample-average trajectories depicted in 
Figure 8. Convergence to stationary values is easily seen. 
Moreover, from the ergodicity property, we can even state that 
these fixed points for the sample-averages are globally attracting 
and that they coincide with the asymptotic time-average (dashed 
lines). The steady-state average values for the proteins A, R and C 
are given by 222.1797, 534.8853 and 549.7195, respectively. 

p53 model 

Let us consider one of the oscillatory p53 models of [67], which 
is described by the reactions 



0 - 

S 3 



Si 

s 2 - 1 



0, Si 
S3, Si 



/(S.,S 3 ) 



(38) 



So. 



where Si is the number of p53 molecules, S2 the number of 
precursor of Mdm2 molecules and S3 the number of molecules of 

Mdm2. The function / (x,y) = ^rjjjr implements a nonlinear 
feedback on the degradation rate of p53. We have the following 
result: 

Result 21 For any positive values of the rate parameters, the oscillatory 
p53 model (38) is ergodic and all the moments are bounded and globally 
converging. 

Lotka-Volterra model 

We consider here the stochastic reaction network 



Stochastic SIR model 

We consider here the following SIR-model, similar to the one in 
[64], defined as 



S + Sj 



Sj, Si Si -f~ Sj 

Si, Si - 0 



(39) 



0 * 


s. 


0 


1 


1, 


I A 


0, 


R 


j'r 


0, 




R, 


R 


* rs 


S. 



S - 0 
S + I k S. 21 (37) 



where birth and death reactions represent people entering and 
leaving the process, respectively. The only bimolecular reaction is 
the contamination reaction which turns one susceptible person 
into an infectious one. The two last reactions represent how 
infectious people are recovering and how recovered people 
become susceptible again. We then have the following result: 

Result 19 For any positive values of the rate parameters, the SIR-model 
(37) if ergodic and all the moments are bounded and globally converging. 

Circadian clock 

Let us consider the circadian oscillator of [65], depicted in 
Figure 6, which is a network involving 9 species and 18 reactions. 
Applying the developed theory on this model, we obtain the 
following result: 

Result 20 For any positive values of the rate parameters, the circadian 
clock model of [65] is ergodic and all the moments are bounded and globally 
converging. 

Using, for instance, the values of [65] and solving for the 
optimization problem (20) using linprog and Yalmip [66], we find 
that ci =402.5768 and cj = 0.1992. Typical trajectories for the 




Figure 5. A -gene repressilator. 

doi:10.1371/journal.pcbi.1003669.g005 
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which is an open analogue of the deterministic Lotka-Volterra 
system of [68]. The first set of reactions represent immigration, the 
second one reproduction, the third one competition due to 
overpopulation and the last one deaths/migrations. We obtain 
then the following result, which is a stochastic analogue of the 
results in [69] obtained in the deterministic setting: 

Theorem 22 Let us define T(v) = [v/y^] and assume that one of the 
following conditions hold: 

1. there exists v>0 such that the matrix T(v) + r(v) T is positive 
definite; 

2. there exists v>0 such that the r(v) + T(v) T is copositive, i.e. 
x r (r(v) + r(v) T ).v>0 for all x>0, and /J,-<5,<0 for aU 
i= 1, . . . ,n. 

Then, the stochastic reaction network (39) is ergodic and all the moments up 
to order [_ 1 + — 2 are bounded and globally converging. 

Schlogl model 

In order to illustrate that the method can be applied to systems 
with more general mass-action kinetics, we consider the stochastic 
version of the well-known Schlogl model [70]: 



2S 



3S 



/o 



k 4 X B 



2S 



(40) 



where S is the main molecule in the network. The above model is 
derived in the supplementary material SI where we have assumed 
that the other molecular populations do not vary over time. Note 
that in the present form the model has an infinite state-space and 
involves a single trimolecular reaction. We then have the following 
result. 




Theorem 23 For any positive values of the rate parameters 
ki,lc2,k3,k4 and any positive values for Xa and Xb, the Markov process 
describing the Schlogl model (40) is exponentially ergodic. 

Note, however, that we cannot say anything on the stability of 
the moments (besides the fact that the first order-moment 
converges) since the condition DD2 does not hold here due to 
the presence of a cubic term. Note that extending the condition 
DD2 to handle more general cases, such as this one, might be 
possible. 

Discussion 

The central theme of this paper is to verify the ergodicity and 
moment boundedness of reaction networks in the stochastic 
setting. Note that even though we mainly consider mass-action 
kinetics in this paper, the framework also applies to more general 
kinetics described, for instance, by Hill functions (see the examples 
on the repressilator and the stochastic switch) and more general 
mass-action kinetics. These results have several interesting and 
important biological implications. 

For example, the ergodicity of a network shows that population- 
level information could be obtained by observing a single 
trajectory for a long time. Such an insight can be used to leverage 
different experimental techniques for a given application. For 
example, consider a clonal cell population with each cell having a 
gene-expression network that is ergodic. Then the stationary 
distribution (at the population level) of the species involved in this 
network can be ascertained by observing a single cell over time. In 
other words, to obtain stationary distributions one can either 
collect samples over time from a single cell (e.g. using time-lapse 
microscopy) or one can take a snapshot of the entire cell 
population at some fixed time (e.g. using flow-cytometry). Due 
to ergodicity, both these approaches will yield the same 
information. Hence, far from being a technical condition, 
ergodicity can have far reaching experimental implications. 

As a property of a network, ergodicity also sheds important light 
on the long range behaviors that can be exhibited by that network. 
One may expect that most endogenous biochemical networks to 
be ergodic in order to achieve robustness with respect to variability 
in initial conditions and kinetic parameters, thus ensuring proper 
biological functions in spite of environmental disturbances. As also 
mentioned in the introduction, ergodicity is a non-trivial property 
which needs to be carefully established and cannot be generically 




Figure 6. Circadian clock model of [65]. 
doi:1 0.1 371 /journal.pcbi.1 003669.g006 



40 50 60 
Time [hours] 



Figure 7. Sample-path of the species of the circadian clock 
model. 

doi:10.1371/journal.pcbi.1003669.g007 
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Figure 8. Time evolution of the sample averages of the species A (top), R (left) and C (right) of the circadian clock model (2000 cells 
averaging). The dashed-lines correspond to the (asymptotic) time-average. 
doi:1 0.1 371 /journal.pcbi.1 003669.g008 



assumed. To illustrate this, let us consider a simplified version of 
the model of carcinogenesis considered in [7 1] which is given by 

h k \2 
0 - Si, Si - S 2 

k 2\ fix) 

S 2 Si, S 2 — 0 

where f(x)= , a>0. When k\ >}>2> the trajectories of the 

cc + x 2 

species grow unbounded, as shown in Figure 9, emphasizing then 
non-ergodicity of the model for this choice of parameters. 

The ideas we use for analysis can also be applied for rationally 
designing circuits in synthetic biology, where it is important 
that the network be (structurally) ergodic in order to ensure that 
the dynamics has the desired behavior irrespective of the initial 
conditions. Such a design is crucial because the initial conditions 
are usually unknown or difficult to control at certain times, 
e.g. after cell division or after the transfection of plasmids in the 
cell. 

Our results on boundedness and convergence of statistical 
moments enable verification of the suitability of a stochastic model 
and to characterize the properties of its steady-state distributions, 
even if such a distribution is not explicitly computable. One 
application of this is to provide justifications and insights for using 
moment closure techniques which have been extensively used to 
study stochastic chemical reaction networks. Some of these 



techniques [72,73] are based on manipulations of the moment 
generating function of the underlying stochastic process. The 
existence of this moment generating function is implicidy assumed 
in such techniques but it may not always hold, thereby 



100 




Time 



Figure 9. State trajectories of the carcinogenesis model (41) 
with the parameters ki—5, = 1, kzi = l, )'2 = 4 and a=l. The 
dashed lines correspond to the average trajectories computed over 
1000 cells. 

doi:1 0.1 371/journal.pcbi.l 003669.g009 
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jeopardizing the validity of the technique. In this article, we show 
that under certain conditions, the distribution of the stochastic 
process is uniformly light-tailed, which proves that the moment 
generating function exists for all time. Certain moment closure 
techniques (see [74,75]) prescribe ways to approximate higher 
order moments as a function of lower order moments. Such an 
approximation is, however, only reasonable if the higher order 
moments are bounded over time. This can be easily assessed with 
our approach and one can even quantify the error by explicitly 
computing the moment bounds as described in this article. 

Finally, the techniques developed here will prove invaluable for 
designing synthetic biological control systems and circuits whose 
objective is to steer the moments of the network of interest to a 
specific steady-state value. Until now, no theory has provided 
guidance for such a design. The specifics are outside the scope of 
this article and will be pursued elsewhere. 
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